/*******************************************************************************
																				
	DESCRIPTION:  	This do file produces graphs that depict the heterogeneity 
					by  duration dependence in our sample.

*******************************************************************************/

clear all
global id_code 119_5
pause on
set seed 2110

*******************************************************************************
* 1. Duration dependence
*******************************************************************************

* Import Duration Dependence Data from 112_4:
use "${data}/119_3_DurationDependenceBetas_Full_2006.dta", clear

* Generate quintiles
xtile quintBeta1 = beta_1, nq(5)
xtile quintBeta0 = beta_0, nq(5)

* Generate empty variables to store the results:
gen t = .

forval i = 0/1 {
	forval j = 1/5 {
		gen f_q`i'_`j' = .
	}
}

gen f_ave = .
gen exp = .

* Loop over which coefficients to use (raw or shrunken):
foreach coef in _shrunk_ind {

	* Calculate E(exp(beta_0)) for the left panel:
	replace exp = exp(beta_0`coef')
	sum exp if !missing(beta_1`coef')
	local mean_exp_beta_0 = r(mean)

	* Calculate E(exp(beta_0)) x E(exp(beta_1 x d) | Q) using quintiles of
	* beta_1 for the LHS panel:
	local i = 1
	
	foreach t of numlist 0(0.2)12 {
			
		replace exp = `mean_exp_beta_0' * exp(beta_1`coef' * `t')
		
		forval q = 1/5 {
			sum exp if quintBeta1 == `q'
			replace f_q1_`q' = r(mean) in `i'
		}
		
		replace t = `t' in `i++'
	}

	* Calculate E(exp(beta_0 + beta_1 x d) | Q) using quintiles of beta_0 for
	* the RHS panel:
	local i = 1

	foreach t of numlist 0(0.2)12 {
			
		replace exp = exp(beta_0`coef' + beta_1`coef' * `t')
		
		forval q = 1/5 {
			sum exp if quintBeta0 == `q'
			replace f_q0_`q' = r(mean) in `i'
		}
		
		replace t = `t' in `i++'
	}

	* Finally, produce overall average profile E(exp(beta_0 + beta_1 x d)):
	local i = 1

	foreach t of numlist 0(0.2)12 {
			
		replace exp = exp(beta_0`coef' + beta_1`coef' * `t')
		sum exp
		replace f_ave = r(mean) in `i'
		replace t = `t' in `i++'
	}

	* Put averages in new frame:
	frame put t f_q* f_ave, into(durDep`coef')
	frame change durDep`coef'
	keep if t != .


	* LHS plot:
	graph twoway (line f_q1_1 t, color(orange_red*0.4) lwidth(medthick)) ///
		(line f_q1_2 t, color(orange_red*0.7) lwidth(medthick)) ///
		(line f_q1_3 t, color(orange_red*1) lwidth(medthick)) ///
		(line f_q1_4 t, color(orange_red*1.3) lwidth(medthick)) ///
		(line f_q1_5 t, color(orange_red*1.6) lwidth(medthick)) ///
		/* (line f_ave t, color(black) lwidth(medthick)) */, ///
		yline(`=f_q1_1[1]', lpattern(dash) lcolor(gray)) ///
		xtitle("Time (months)") xscale(range(0 12)) xlabel(0(3)12) xtick(0(1)12) ///
		legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") position(6) rows(1) ///
			symxsize(*0.5) size(medsmall) linegap(2) ///
			region(margin(small)) subtitle(" ", size(minuscule)) ///
			title("Quintiles of {&beta}{sub:D}, Fixing {&beta}{sub:0}", ///
				size(medsmall) alignment(top) color(black))) ///
		ytitle("") ylabel(0.2(0.1)1) ///
		/// t2title("{bf: E}[{&pi}{sub:d}|{bf: E}({&beta}{sub:0}), Q{sub:D}]") ///
		/// t2title("A. By Quintiles of Duration Dependence, Fixing Initial JF") ///
		graphregion(color(white)) ///
		name(g1, replace)

	* RHS plot:
	graph twoway (line f_q0_1 t, color(ebblue*0.4) lwidth(medthick)) ///
		(line f_q0_2 t, color(ebblue*0.7) lwidth(medthick)) ///
		(line f_q0_3 t, color(ebblue*1) lwidth(medthick)) ///
		(line f_q0_4 t, color(ebblue*1.3) lwidth(medthick)) ///
		(line f_q0_5 t, color(ebblue*1.6) lwidth(medthick)) ///
		/* (line f_ave t, color(black) lwidth(medthick)) */, ///
		yline(`=f_q1_1[1]', lpattern(dash) lcolor(gray)) ///
		xtitle("Time (months)") xscale(range(0 12)) xlabel(0(3)12) xtick(0(1)12) ///
		legend(order(1 "Q1" 2 "Q2" 3 "Q3" 4 "Q4" 5 "Q5") position(6) rows(1) ///
			symxsize(*0.5) size(medsmall) linegap(2) ///
			region(margin(small)) subtitle(" ", size(minuscule)) ///
			title("Quintiles of {&beta}{sub:0}", ///
				size(medsmall) alignment(top) color(black) ///
				) ///
			) ///
		ytitle("") ylabel(0.2(0.1)1) ///
		/// t2title("{bf: E}({&pi}{sub:d}|Q{sub:0})") ///
		/// t2title("B. By Quintiles of Initial JF") ///
		graphregion(color(white)) ///
		name(g2, replace)
		
	* Combine graphs:
	graph combine g1 g2, rows(1) ycommon ///
		l1title("Average 6-Month JFR", size(medium)) ///
		graphregion(color(white)) ysize(4.5) xsize(8) name(g3`coef', replace)
	
	* Store them:
	if "`coef'" == "_shrunk_pop" {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_DurationDependence_ShrunkPop.pdf", as(pdf) replace 
	}
	else if "`coef'" == "_shrunk_ind" {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_DurationDependence_ShrunkInd.pdf", as(pdf) replace 
	}
	else {
		graph export "${output}/${id_code}_Visualizing_Heterogeneity_DurationDependence.pdf", as(pdf) replace 
	}
	

	* Go back to default frame:
	frame change default
}